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Abstract 

A novel adaptive technique for electromagnetic Particle In Cell (PIC) plasma 
simulations is presented here. Two main issues are identified in designing adap- 
tive techniques for PIC simulation: first, the choice of the size of the particle 
shape function in progressively refined grids, with the need to avoid the exertion 
of self-forces on particles, and, second, the necessity to comply with the strict 
stability constraints of the explicit PIC algorithm. The adaptive implementa- 
tion presented responds to these demands with the introduction of a Multi Level 
Multi Domain (MLMD) system (where a cloud of self-similar domains is fully 
simulated with both fields and particles) and the use of an Implicit Moment 
PIC method as baseline algorithm for the adaptive evolution. Information is 
exchanged between the levels with the projection of the field information from 
the refined to the coarser levels and the interpolation of the boundary condi- 
tions for the refined levels from the coarser level fields. Particles are bound to 
their level of origin and are prevented from transitioning to coarser levels, but 
are repopulated at the refined grid boundaries with a splitting technique. The 
presented algorithm is tested against a series of simulation challenges. 

Keywords: Particle-In-Cell, Adaptive, Implicit, Particle Splitting 



1. Introduction 

The fascinating task of plasma physics simulations is cursed with the co- 
existence of multiple space and time scales in all relevant phenomena. As an 
example, Fig. [l] reports the typical scales of the plasmas populating different 
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regions of the Earth magnetotail. Notice that the smahest scales are of the 
order of the tens of microseconds and of the hundreds of meters, considerably 
smaller than those of interest for the most important processes observed in 
space, such as the geomagnetic storms caused by the arrival of Coronal Mass 
Ejections (CME) at Earth or the substorms originated in the Earth environment 
by magnetic reconnection in the magnetotail region|lJ. Those events involve the 
whole Earth environment and last for minutes to hours to days, making their 
simulation impossible if all plasma scales are to be encompassed. 
A common strategy to make these simulations feasible with the computational 
resources currently available is to neglect the smaller, kinetic scales and use 
fluid[2] of hybrid (fluid electrons, kinetic ions) [21 [3] approaches, at the expense 
of complete consistency and of the kinetic effects. 
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Figure 1: Representation of the typical plasma scales observed in the Earth 
magnetotail. Notice the wide range of spatial and temporal scales which makes 
the simulation of the relevant phenomena happening in such a system an ex- 
tremely expensive task under the computational point of view. 



If such effects need to be retained, explicit Particle In Cell (PIC) meth- 
ods ^ are one of the most frequently employed options. However, their strict 
stability constraints (Ref. [1] and [S]) impose notable consequences for space sim- 
ulations [B], requiring for the smallest Debye length (an intrinsic plasma scale 
measuring the shielding length in a plasma |SJ|1]) and the fastest electron plasma 
frequency to be resolved. This task is particularly challenging since often space 
plasmas present regions of high density in contact with regions of lower density, 
resulting in a wide range of Debye lengths: for such a simulation not to incur 
into stability issues, the grid spacing needs to be chosen to resolve the smallest 
Debye length in the system. Adaptive schemes can be devised to resolve in each 
region only the local Debye length and the local plasma frequency, resulting 
in considerable savings [7j: the use of adaptive grids resolving the local Debye 
length allows to save computational resources in the areas where high resolu- 
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tions are not required. The development of adaptive algorithms in the context 
of PIC methods is however complicated by the fact that PIC methods are es- 
sentially based on the assumption that the particles have a finite size which is 



fixed and equal to the cell size (see later Sec. 2.1), thus making the approach 
naturally suited to uniform grids and of non-straightforward application to non 
uniform grids. Nevertheless there are noteworthy attempts to expand the PIC 
method to non-uniform grids, as reviewed below for some of the most notable 
approaches of interest to space weather problems. 

We organize such a review in three groups. First, we consider Moving Mesh 
Adaptation (MMA) methods, where only one grid level is simulated (usually, 
adaptive grid methods prescribe the simulation of multiple grid levels), with a 
non-uniform structure and a connectivity that allows one to map the simulated 
grid to a uniform logical grid. The non uniform grid is regenerated at each time 
step as the end result of a process that distorts it by attracting more points in 
the regions of interest and away from regions where the system is smooth and 
uneventful. This approach has been followed for example within the Celeste 
Implicit Moment method and is reviewed in Ref . [HI [HI HH] ■ 
Second, Adaptive Mesh Refinement (AMR) [TT] can be used for the field part 
of the PIC method. Notable examples of this approach, which was developed 
especially for space weather application and has been shown to resolve the lo- 
calized regions of reconnection embedded in larger systems [7] , are Ref. [TH [7] . 
Lastly, unstructured grids [T3] can be used starting from finite element methods 
to discretize the fields. 

All these approaches need to deal with the fact that, as the grid spacing changes 
locally, the particles and the cells cannot always have the same size. The PIC 
derivation is based on such an equality to provide a simple efficient interpola- 



tion between particles and cells (see, for example, Ref. [13] and Sec. 2.1 later 
in the paper). The first two approaches relax the requirement that the particle 
size remains constant and base the interpolation on the local grid spacing: the 
particles have the local size of the cells they are embedded in. For the MMA 
approach this is achieved by assuming that the particle shape function is con- 
stant in the logical space ^ rather than in the physical space x. If the shape 
functions are chosen in the logical space, the interpolation functions can still be 
computed with the b-spline chain rule and the same formulas of uniform grids 
can still be used. In the second group of adaptive techniques, a patch-based 
AMR approach is typically used: the particles have the size of the cells of the 
patch where they reside and interpolation can proceed in each patch as if it was 
a uniform grid. 

This simplification of the interpolation step on adaptive PIC has a serious con- 
sequence: the exact derivation of Ref. [6] breaks down and the changing particle 
shape introduces new terms proportional to the temporal derivative of the shape 
function. The derivation of the PIC method outlined in Ref. [6 and briefiy re- 
minded later in the paper assumes a fixed particle shape: a simple analysis 
proves that, within this assumption, momentum is conserved exactly when the 
same interpolation function is used for all quantities [3 [15] and provided that 
the solution of the Maxwell's equations also conserve momentum (i.e. the nu- 
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merical Green's function reflects the properties of anti-symmetry of the exact 
continuum Green's function). When the shape functions for particles change in 
time, instead, additional terms related to the temporal variation of the shape 
function have to be added to the equations of motion. However, in practice 
these additional terms are neglected and the resulting schemes do not conserve 
momentum, with the additional side effect that particles experience a self-force 
due to their traversing regions of different sizes and thus to the change in time 
of their shape function [TBI Hi] • In the MMA approach, this error is controlled 
by using smooth grids where the spacing changes progressively: since the ne- 
glected terms are proportional to the changes in the shape function, when the 
changes are small the error is small. However, in the AMR approach the change 
in particle shape is sudden at the interface between patches and leads, first, to 
the necessity of refining the grids in small steps (refinement factors of two are 
the usual choice) and, second, to severe errors that need to be corrected with 
appropriately devised methods [T71 [TB] . 

In this paper, a novel approach which allows to avoid such inconveniences is 
presented. 

The following two points are the main elements of discontinuity with respect to 
the existing implementations: first, while in the previous works the grid area 
needing an increased resolution was simulated completely with fields and parti- 
cles only at the most refined level, here also the so-called "coarser grids" (i.e., 
the grids with larger spacial resolutions) are retained and fully simulated along- 
side the more refined areas, thus making this implementation a Multi Level 
Multi Domain (MLMD) approach. All the simulated domains are self-similar 
and fully functional, each of them being evolved in time for both fields and par- 
ticles, which are bound to the birth level and whose shape function is tailored 
on the local grid size. Second, the baseline algorithm used here is the Implicit 
Moment PIC Method (see Sec. 3.2), in contrast with the explicit PIC used in 
the other cases. 

The structure of the paper intends to guide the reader towards the motivations 
that lead to such an implementation. 

First, in Sec. [2j the challenges to face in the development of adaptive PIC 
methods are recalled and pinpointed in the necessity to adapt the method to 
the variation of the grid size (Sec. 2.1 ) and in the stability dictated limits on the 
grid resolution required by explicit PIC methods (Sec. 2.2 1. In Sec. [s] our solu- 
tions to such challenges are presented: using multiple, complete and self-similar 
domains (i.e., a MLMD approach) offers an elegant solution to the problem of 
adapting the particle shape functions to the changing grid resolution (Sec. 3.1), 
while the Implicit Moment PIC algorithm allows to notably relax the above- 
mentioned stability conditions (Sec. 3.2). In Sec. |4] the approach proposed is 
presented in its conceptual guidelines, while in Sec. [5] the details of the infor- 
mation exchange between the levels are provided. Sec. |6] reports the results of a 
battery of tests performed in a ID environment: a set of challenges is selected 
to address specific issues arising in a MLMD implementation and provides the 
occasion for a series of comments on fields and particles coupling in the system. 
Finally, conclusions are drawn in Sec. [7] 
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2. Challenges for adaptive Particle In Cell methods 



In this Section the two main challenges to be faced when developing adaptive 
PIC methods are recalled. However, to fully appreciate them a reminder of the 
derivation of the PIC method is needed. 

The PIC method is a mathematically rigorous discretization in a Lagrangian 
form of the Vlasov-Maxwell system describing a plasma at the kinetic level [5] . 
The phase space distribution function /^(x, v,t) for a given species s (electrons 
or ions) is defined as the number density per unit element of the phase space 
and is governed by the Vlasov equation: 

^+v.f^ + ^(E + vxB).f^.O (1) 
at ox TTis ov 

where qs and rus are the charge and mass of the species, respectively. The 
derivatives are written in the standard vector notation in the three-dimensional 
configuration space x and three-dimensional velocity space v, forming together 
the classical phase space. 

The electric and magnetic fields are given by the Maxwell's equations, 

V X E = — — 
dt 

eoV • E = p 
V - B = 0, 

where the sources of the Maxwell's equations are the first two moments of the 
particle distribution function, the charge density, 

p(x,t) = Vg, / /,(x,v,t)dv (3) 
Jy 

and the current density, 

j(x,i) = Vg, / v/,(x,v,t)dv (4) 
Jy 

The integrations that define the moments are done over the velocity space V. 
The two key founding points of the PIC method that prevent its straightforward 
application in an adaptive grid framework are the following. 
First, the size of the computational particles of the classical PIC method is taken 
to be constant in time and equal to the grid size. Transitioning to an adaptive 
grid would seem to require to sacrifice one of the two: either the particles have 
a constant shape that no longer fits the grid size as the particle moves in regions 
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of different resolution or the particles size is always the local grid size, at the 
sacrifice of the constancy of the particle size. The constancy of the particle size 
is a cornerstone of the derivation and, if sacrificed, raises fundamental issues for 
particle motions, while the equality of particle and grid size is of great practi- 
cal convenience in the implementation of the PIC method: giving it up means 
a much more cumbersome algorithm for grid-particle projection and interpola- 
tion. The approach presented here does not give up on cither aspect and retains 
constant particle sizes while maintaining the size of particles and cells equal. 
Second, PIC methods are strongly bound by stability limits: the grid spacing 
cannot exceed the local Dcbyc length (with a factor of order 1) regardless of 
the scales of interest, thus significantly limiting grid adaptivity. This seemingly 
reduces adaptive PIC to methods to adapt the grid to the local Debye length, 
still a very \iscful tool in systems where the Debye length changes over a sig- 
nificant range, as it is indeed the case in many space and laboratory plasmas. 
However, this limitation should be removed for a truly adaptive PIC where the 
local grid spacing is chosen solely based on the accuracy of the resolution de- 
sired: consider that in many systems, in space and laboratory, quasineutrality is 
very strongly satisfied and nothing of interest happens at the Debye length. The 
method proposed bypasses the stability constraints of the explicit PIC method 
by adopting an Implicit Moment PIC method as the baseline algorithm to evolve 
in an adaptive direction. In the Sections below, the two apparent roadblocks 
mentioned above and the approaches proposed to overcome them are discussed 
in more details. 

2.1. Changing grid and particle size 

In PIC methods, the phase space is discretized into a collection of super 
or computational particles, each representing a number of physical particles 
close to each other in phase space. As a consequence of this assumption, the 
distribution function of each species s is written as: 

/s(x,v,t) = ^/p(x,v,i) (5) 

p 

where the index p spans the computational particles of species s. The specific 
distribution function fp assigned to each computational particle has a number 
of free parameters which have the physical meaning of position and velocity of 
the computational particle and whose time evolution determines the numeri- 
cal solution of the Vlasov equation. The distribution function fp{x,v,t) of a 
computational particle is thus assumed to be the following: 

fp{K,v,t) = 7Vp5x(x- Xj,(i))^v(v - Vp{t)) (6) 

where and 5v are the shape functions in space and velocity of the compu- 
tational particle and Np is the number of physical particles that are present in 
that element of phase space. The common choices for the velocity and spatial 
shape functions are respectively Dirac deltas 
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S'v(v - Vp) = S{vi - Vi^p)5{v2 - V2,p)S{v3 - V3^p), (7) 

and b-splines [18^, usually of order (Cloud In Cell, CIC methods), 



5x(x - xp) = / , 6, ( V'^'O ( f . 

(8) 

where Axi^p, Ax2,p and Ax^^p are the lengths of the support of the compu- 
tational particles in each spatial dimension. 

When integrating the shape functions over the grid, i.e., when performing the 
following integral in all directions 

Sx{x-Xp)d-x = J Sx{x - Xp)bo (^ ^^^^ ^ dx, (9) 

the interpolation functions W{'x.g — Xp) play a fundamental role. They are 
in fact used both to calculate particle moments on the grid points, as in 

^9 p 

(10) 

jf = T7 51 '?pVpTy(Xg - Xp), 

^9 p 

where Vg is the cell volume, and to calculate the fields acting on particles 
from the fields defined on the grid, as in 

Ep = ^Egl^(Xg-Xp) 

9 

(11) 

Bp = 5lB3W(xg-Xp). 

9 

In the cases when the support for the spatial shape function is equal to the 
grid spacing, Axi^p = Axi^g, and the shape function is a b-spline of order 0, the 
interpolation functions simply reduce to 



Wi.. - X,) . (^^) Oh, (a^) (a^) . (12) 

As already mentioned, issues arise if Ax^^p 7^ Axi^g (one of the possible 
solutions to the issue of particle shape in an adaptive grid, with particles of fixed 
shape moving across grids with different resolutions), since the mathematical 
description for the interpolation function becomes more complicated. 
If the shape function does not explicitly include a time dependence, that is if 
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= 0, 



(13) 



dt 



the equations for the motion of the computational particles in the PIC 
method follow directly from the Vlasov equation by simply noting that 



where Xi represents the three space variables. They are thus identical to those 
of a particle with center of mass properties given by the parameters of the 
computational particle, Xp and Vp. If the terms dSxJdt are present, instead, 
as in the case of particle shapes being tailored on changing grid sizes (the most 
used solution in AMR PIC), the equations also include additional terms that 
make them not resemble anymore those of physical particles. In particular, they 
would include self-forces [121 EZ] that act on the superparticles in consequence 
of the change in the particle shape. 

One possibility to deal with the issue would be to limit the value of the terms 
dSxi I dt with respect to the other terms of the equation of motion: if the particle 
shape changes slowly with the particle motion, the particle shape remains, in 
first approximation, an adiabatic invariant of motion. The equations for the 
superparticles remain then almost identical to those of real particles and the 
self forces become negligible. This approach, which requires the grid resolution 
to change slowly in space, has been followed in the variational grid adaptation 
[HI [19] and has been implemented recently in the code DEMOCRITUS [TO] . 
The need to limit dS^Jdi is incidentally one of the reasons why AMR PIC 
methods tend to prefer small resolution jumps between coarser and refined grid, 
the typical value being two. 

Another possibility is to split particles that move from regions with lower 
resolution to regions with higher resolution, and, similarly, coalesce particles 
traveling in the opposite direction. Reliable methods have been developed to 
coalesce and split particles in order to force the particle population to maintain 
its total number per cell constant and have been successfully applied, for exam- 
ple, in Ref. [7|. The method does not attempt to exactly remove the self- forces 
created by the violation of the condition in Eq. |13[ but rather eliminates them 
in a statistical, average way, compatibly with the statistical nature of the PIC 
approach. Note however that, while the splitting procedure for particles is rela- 
tively straightforward, as described in Ref. [50], the inverse coalescence process 
is less immediate to code (particles with the same mass and charge and close 
in phase space have to be identified) and does not exactly conserve the total 
energy and the local distribution function [^. 

Finally, the exact removal of the self-forces can be achieved by giving up on 
their constant shape and complicating the interpolations and grid projection al- 
gorithms. This approach has been proposed within unstructured meshes meth- 
ods [T3] by simply computing the complex geometric overlap of a constant parti- 
cle shape with the local grid arrangement. In AMR methods the grid is uniform 



dxi. 



(14) 



dt dxi 



dt 
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in patches and the particle shape changes only at the boundary between re- 
gions of different resolution. There, the condition of Eq. [13] is violated and a 
correction is applied. 

2.2. Stability dictated limits on grid resolution 

The stability dictated limits for explicit PIC methods and their consequences 
for plasma simulations have already been commented in Sec. [T] and [2] Let us 
recall here the origin of the necessity to resolve the fastest electron plasma 
frequency and the smallest Debye length in a system simulated with an explicit 
PIC method. 

When an explicit formulation of the coupling between particles and fields is used, 
the particles are moved at each time step in the fields known at the beginning of 
the time step and the fields evolve for the same time step based on the recently 
updated and fixed particles. This works only if the time step is small, thus 
introducing the stability condition 

UJpeAt < 2. (15) 

The other stability condition, cAt < Ax, stems from the explicit discretization 
of the Maxwell's equations and is normally not the main concern in space plas- 
mas, since it can be easily removed with semi-implicit methods that discretize 
implicitly only the Maxwell's equations; the most challenging part is to keep 
the particle- field coupling implicit. 

A more serious drawback from the perspective of adaptive methods is that 
explicit PIC algorithms are subject to the finite grid instability. The mathemat- 
ical study of the finite grid instability requires an analysis of all computational 
steps with the Laplace transformation of time and the Fourier transformation 
of space (see Ref. 0] E] and citations therein) ; the summary of the analysis is 
that to avoid the finite grid instability in explicit particles methods, the grid 
spacing must be chosen to satisfy the constraint 

Ax/Xoe < <; (16) 

where A De is the Debye length and the parameter is of order one and depends 
on the details of the implementation. For the CIC scheme, the literature reports 

The practical consequence of this numerical instability is a tremendous nu- 
merical heating of the plasma characterized by an alternatively positive-negative 
variation of the electric field accompanied by a correlated zig-zag perturbation 
of the phase space. Within a few cycles, the energy reaches the bounds of over- 
flow in the machine representation of numbers. The finite grid instability must 
thus be avoided at all costs, requiring its stability constraint (ujpeAt < 2 and 
Ax/Xue < ^) to be respected everywhere in the domain for all directions for 
the fastest Up,, and the smallest Xoe in the system. Indeed, a practice advice is 
to use a considerably smaller time step of order uipeAt — 0.1 and Aa; < \De^ to 
avoid numerical heating 4J. 



9 



It has already been observed that an AMR approach helps reducing the 
impact of the explicit PIC stability constraints in cases where the processes of 
interest are at the Debye length level but the presence of regions of higher density 
in contact with regions of lower density results in a large range of Debye lengths. 
However, less fortunate situations may reduce the benefits granted even by the 
application of a AMR algorithm to an explicit PIC. Consider, for example, that 
most reconnection events in space are dominated by processes that are on the 
electron inertial length, de = c/wpe. This scale is still two orders of magnitude 
larger than the Debye length: having to resolve it purely for numerical stability 
reasons implies a waste of two orders of magnitude in each direction, a million 
times more cells than necessary [51 [H] . 

If using uniform grids, this constraint thus requires to use cell sizes several 
orders of magnitude smaller than the scales of interest. However, also if adaptive 
grids are used the consequences of the finite grid instability on explicit PIC 
simulations are heavy, since the range within which the refinement factor of the 
refined grids can be chosen is severely limited. Indeed, the coarsest cells still 
need to resolve the local Debye length even in regions where there is no interest in 
having resolutions that small. Taking again the example of the regions in Fig. [l] 
to study a reconnection event in the magnetotail the ideal situation would be to 
resolve up to the electron skin depth in the vicinity of the reconnection region, 
but up to the ion skin depth or to even larger scales outside of it. Unfortunately, 
that is not possible even using adaptive explicit PIC. 



3. Tackling the challenges 

In this Section, the challenges to the development of PIC methods for adap- 
tive grids are tackled and solutions are proposed. A MLMD approach is used to 
bypass the issue of the size of the particle shape functions in adaptive grids and 
results into introducing additional benefits to the system, such as the possibility 



of sudden jumps in the Refinement Factor and ease of programming (Sec. 3.1). 
The strict stability conditions for the explicit PIC methods, so harmful also in 
the prospect of an adaptive grid evolution, are instead dealt with by switching 



to an Implicit Moment PIC approach (Sec. 3.2) 



3.1. Multi Level Multi Domain approach 

The adaptive approaches conceptually closest to the present implementa- 
tion [ini m] choose, for the sake of performances, to simulate computational 
particles only at the most refined levels and to assign them a shape function 
with dimensions equal to the most refined grid spacing. This solution raises the 
issue of self-forces and the need for corrective actions jl7) . 
In this paper, a different, MLMD solution to the particle size dilemma is pro- 
posed: instead of relegating particles to the most refined levels, the system is 
simulated as a cloud of overlapping domains, each logically complete in both 
fields and particles and fully identical to the others, except in position, resolu- 
tion and size. 
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Particles are created at initialization in each level with the size of the spatial 
shape function of Eq. [8] equal to the spacing of the grid they belong to and are 
bounded to their level of origin, not being allowed to transition from coarser 
to refined grids and vice versa. The details about how particles are dealt with 
at the boundaries between the coarser and refined grids are provided later in 
Sec. |5.3[ for the purpose of the present discussion it is sufficient to mention that 
refined grid particles are lost when they exit their domain of origin and new 
refined level particles are created at the boundary cells of the refined grids with 
the splitting algorithm devised in Ref. [20] and already used in Ref. [21]. Such 
an algorithm replaces a coarser grid particle with a collection of refined grid 
particles whose shape functions collectively reproduce the spatial extension of 
the original space function, thus eliminating self forces at the boundary in a 
statistical way. 

This apparently simple MLMD solution grants a number of advantages. The 
fact that the shape function of particles has the same size of the grid granularity 
means that the simple interpolation functions of Eq. [12] can be used to straight- 
forwardly calculate the sources of Maxwell's equations (Eq. [lO] ) and the fields 
acting on a particle (Eq. [11[ ) : no complicated interpolation techniques need to 
be devised. More importantly, the particle shape function staying constant in 
time implies that no self forces are exerted over the particles, thus eliminating 
the need to correct for them. Lastly, since a refined area is simulated also with 
the coarser resolution both in fields and particles, no coarser grid particles have 
to be created by coalescence of the refined grid ones, thus avoiding the prob- 
lems of conservation of total energy and local distribution function reported in 
Ref. [7]: only the "safest" part of the algorithm of Ref. [SCT is used here. 

A graphical depiction of the proposed system is given in Fig. [2] an arbitrary 
number of levels with different Refinement Factors is fully simulated in fields 
and particles. Particles are present at each level with shape function sized on the 
local grid size and are used to calculate locally the sources of the Maxwell's equa- 
tions, which are solved for on each level. Each level completes the same 
operations as the others and the communication between the grids is limited to 
the " downwards" (interpolation of the boundary conditions for the refined levels 
from the coarser level fields) and "upwards" (projection of the refined fields from 
the refined to the coarser levels) operations marked respectively with number 
(1) and (2) in the picture and described in more details in Sec.js] 

Such a system, apart for providing an elegant solution to the shape func- 
tion problem, presents a few more aspects which made the MLMD approach 
attractive to the authors. 

Unlike other AMR approaches, in fact, it removes the need to proceed in 
soft jumps of resolution differential. Typically AMR methods (see, for example, 
Ref. [7]) progressively create cells smaller by a factor of two in each direction, 
due to the increase of the self-forces with increasing difference between the 
particle shape function on the different levels and to the explicit PIC stabil- 
ity constraints. Here (also thanks to the implicit moment implementation, see 
Sec. 3.2 1, such a necessity is eliminated. This is a relevant achievement consider- 
ing that kinetic problems tend to jump in big steps. Consider, for example, the 
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Figure 2: Proposed Multi Level Multi Domain approach. The system is simu- 
lated as a cloud of self-similar, partly overlapping domains all complete in fields 
and particles. Particles are present at all levels, including the overlapping ar- 
eas, and are sized according to the local grid size. Each level executes the same 
operations for fields and particles and the exchange of information for the fields 
is limited to the "downwards" interpolation of the boundary conditions for the 
refined levels from the coarser level fields (1) and the "upwards" projection of 
the refined fields from the refined to the coarser levels (2). Particles are regen- 
erated at the boundaries of the refined grid with a splitting algorithm and are 
lost when they exit the refined areas, thus eliminating the need of coalescence 
operations. 

now familiar magnetic reconnection scenario, with ions and electrons becoming 
unmagnetized at different distances from the X-line: in an ideal world, an adap- 
tive method should allow a Refinement Factor proportional to the square root 
of the mass ratio in transitioning from the ion to the electron diffusion region. 
Moreover, the concept presented is intrinsically simple and object oriented in 
nature and as such the coding of the algorithms is by far easier than most AMR 
codes based on trees or patches. The promise of simplicity is not just in the con- 
ceptual ingredients but also in the operations, which are identically performed 
on all levels, regardless of the level of grid refinement. 

3.2. The Implicit Moment Particle In Cell method 

Implicit methods have been considered for several decades [2H [231 HH 123 
an answer to the problem of the large number of time steps and the high reso- 
lution needed to run realistic problems with explicit particle methods. Here, an 
Implicit Moment PIC approach is adopted in the implementation typical of the 
family of codes started by Venus [53] and continued by Celeste [27] , Parsek [5S] 
and iPicSD [55^. The aim is to bypass the stability constraints of the explicit 
methods and allow more flexibility in the choice of the Refinement Factors. 
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In the presented implementation, the particle mover is based on the 9 scheme [271 
ED]: 

n+l , QsAt ' — 



V = V 

V V 



where a decentering parameter 0, which is used to compute the weighted aver- 
ages between the old and new time level quantities, as in ^1/"+^ = v[/"(i — (?) -|- 
Vf+^^j is used to vary the properties of the scheme. 

The velocity equation is more conveniently rewritten as: 



v;'+i/2 ^ + /33E;'+«(x;'+i/2) (18) 



with Pis = qpAt/2T 



and a" defined as 



(19) 



1 + {PsB- 



{I - f3J X B" + /32b"B") (20) 



and representing a scaling and rotation of the velocity vector. 

Eq. [it] and Eq. [18] for particle motion are in implicit form. In the Implicit 
Moment method, the Maxwell's equations, discretized as in 

b;'+i - b;; = -Aiv X e;'+^ 



(21) 



V • B^'+i = 0, 



are also solved implicitly, introducing a coupling between the equations for par- 
ticle motion and those for the fields. This dependence is removed in the Implicit 
Moment method as summarized in Fig.jsj modified from Ref. |14j . and described 
in Refs. [231 1571 [TH |^ : the sources of the field equations are approximated 
using the moment equations instead of being calculated directly from the par- 
ticles. 

Once the field equations are solved within this approximation, the rest of 
the steps can then be completed sequentially: with the new fields, the particle 
equations of motion can be solved (either with a predictor corrector PC [30] or 
a Newton approach [32]) and the new current and density can be computed for 
the next computational cycle. 
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Moment Eqs. 
J(Eg, Bg) 
N(Eg, Bg) 



Field Solver 
Eg, Bg 



Grid -> Particle 
Eg, Bg -> Fp 



Particle Mover 
Fp->Xp 



Particle -> Grid 
Xp->Ng, Jg 



Figure 3: Blocks needed in a cycle of the implicit moment PIC scheme. Mod- 
ified from Ref. [14j. 



The key step is thus to derive a suitable set of moment equations that can 
approximate the response of the plasma particles to the fields over a compu- 
tational cycle. This is done through a series expansion of the interpolation 
function appearing in the expression for the field sources. The approximation 
is performed with respect to the particle position, choosing as the center of the 
expansion the particle positions at the beginning of the computational cycle (in 
direct implicit methods, instead, the expansion is centered around a guess of 
the new advanced position [25 ): 

W(x-x^+i) = W^(x-xp+(x-xp.VM^(x-xp + ^(x-xp(x-x;) : VVM/(x-xp+. . . 

(22) 

where a tensor notation is used. 

The expansion of Eq. [22] can be used to compute the field sources directly 
using particle information only from the previous computational cycle and re- 
moving the need to iterate over the particle and field equations. The details 
of the simple but tedious algebraic manipulations are provided in Ref. [57] , the 
final answer being: 

where the following expressions were defined: 

is 

Us 



(23) 



= Ep9pVpW^(x-xp 
Ep'7pVpVpH/(x-xJJ) 



(24) 



with the meaning, respectively, of current and pressure tensor based on the 
transformed hatted velocities. The fi^ term is an effective dielectric tensor 
which expresses the feedback of the electric field on the plasma current and 
density: 

m: = -— <■ (25) 
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The expression in Eq. 23 for the sources of the Maxwell's equations provide 



a closure for the Maxwell's equations themselves. When Eq. 23 is inserted in 
Eq. [21] the Maxwell's equations can be solved without further coupling with the 
particle equations, as in the consistent second order formulation from Ref. |33j : 

{cdAtf [-V2E"+'' - VV • (/x" • E"+^)] + e" • E"+^ 

= E" + {c9At) (w X B" - ^j"^ - {c0Atf V47rp", ^^^^ 

where /x" = J^s /^s e" = I + /x", I being the identity tensor. 

As shown in Ref. |33] , the second order formulation for the electric field needs 
to be coupled with a divergence cleaning step to ensure that Gauss's law (Eq. [2]- 



3, discretized in Eq. 21 -3) holds true at each time step. The magnetic field is 



then computed directly once the new electric field is known, as in Eq. 21 1. The 
discretized equations and their boundary conditions form a non-symmetric lin- 
ear system that is solved using the Generalized Minimal Residual (GMRES) [31] 
method. For the divergence cleaning equation. Conjugate Gradient (CG) meth- 
ods or Fast Fourier Transform (EFT) methods on uniform grids j34| can be used 
since the discretized equation leads to a symmetric matrix. 

The above summarized Implicit Moment formulation has the notable advan- 
tage of removing the strict stability limits of the explicit PIC. 
The stability properties of the implicit moment method have been studied ex- 
tensively in the past [21] : the implicit particle mover removes the need to resolve 
the electron plasma frequency and the implicit formulation of the field equations 
removes the need to resolve the speed of light. The time step constraints are 
thus replaced by an accuracy limit arising from the derivation of the moment 
equations using the series expansion in Eq. |22[ This limit restricts the mean 
particle motion to one grid cell per time step ^24) . i.e. 

Vth,e^t/Ax < 1, (27) 

The finite grid instability limit for the explicit method. Ax < <;XDe, is replaced 

by m 

Ax < e-^vth,eAt, (28) 

that allows large cell sizes to be used when large time steps are taken. 

This is the key feature of interest of the implicit method: time and space 
scales can be chosen freely according to the desired accuracy (refer to Refs. [H 
[551 m] for details on how the non-resolved scales are dealt with) , but their ratio 
is not free, since it must stay within the bounds just outlined: 

e < vtK^At/Ax < 1, (29) 

The upper limit is generally the main concern and should never be violated. 
The lower limit, e, is in practical cases usually a very small number that can 
often be approximated by zero. 

The gain afforded by the relaxation of the stability limits is already notable 
in non-adaptive cases, but becomes particularly precious in the adaptive case. 
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As regards the non adaptive cases, the new stabihty constraints remove the 
necessity to resolve the faster electron plasma frequency and the smaller Debye 
length while retaining the sub-At scales and thus the possibility to exchange en- 
ergy between sub-At fluctuations and particles (notice that in other approaches, 
e.g. the gyrokinetic or hybrid one such processes are completely removed 
and the energy channel towards them is interrupted). Instead, when adaptive 
schemes are considered the increased range of time steps and grid spacings avail- 
able means that the level of refinement between coarser and refined levels may 
be chosen more freely, without the concern of resolving the Debye length in all 
grids. As regards the choice of the time steps, two possibilities are now open: 
the time step can either be adjusted to the grid spatial resolution, with the 
refined levels using a fraction of the time step of the coarser levels in order to 
maintain the same Aa;/At ratio across the levels, or At can be set for the entire 
system to the value required by the more refined level, the critical part of the 
inequality of Eq. [29j Here this second solution is preferred, since in a future 
parallel implementation with exchange of information between the levels adopt- 
ing a larger time step for the coarser domains in absence of a sophisticated work 
sharing algorithm would just make the coarser grids processors remain idle while 
waiting for communication from the refined grids processors. Notice again that 
the e < Vth,e^t/Ax part of the stability constraints is, according to practical 
experience, not a significant hindrance and that the value of e can be rather 
small, thus not significantly limiting the spacial grid resolution of the coarser 
grids when the At is chosen according to the refined grids requirements. 
Observe additionally that Implicit Moment methods come with an embedded 
numerically-enhanced Landau damping at high wavenumbers. This characteris- 
tic proves to be a precious ally against the issue of the reflection of low-A waves, 
supported by the refined grids but not by the coarser grids, at the boundaries be- 
tween the levels |16| , a problem which becomes dramatic for explicit algorithms 
if certain stencils are used for the discretization of the fields |36J. 

4. Overview of the Multi Level Multi Domain method 

To summarize the argumentations developed above, the MLMD approach 
proposed here is based on the following principles: 

• the physical system is represented as a cloud of overlapping domains or- 
ganized in levels, each one of them conceptually identical to the others 
except for a scale factor and a translation. At the coarsest (upper level 
in Fig. [2]), the whole domain is discretized in a coarse uniform grid. At 
the next levels, the areas of interest are further discretized by other do- 
mains identical to the upper levels in everything (including number of 
cells and particles) except for being rescaled in size down to a fraction, 
not necessarily in soft jumps; 

• while the top coarsest level covers the whole computational domain, the 
lower levels with higher resolution cover only progressively smaller parts 
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of the physical system, the areas where the enhanced resolution is really 
needed; 

• the scaling affects cells and particles alike, with particles populating all 
the levels, not just the most refined, and having the same size as the cells 
they move in. As a consequence, the particles of the refined levels have a 
smaller shape function than those of the coarser levels; 

• since each level has its own cells and particles and performs the same op- 
erations as the others, from a software point of view all levels are identical 
and in a object oriented approach they are different instantiations of the 
same class. For each level, the instantiation is rescaled and shifted; 

• the communication between the levels consists of the projection of the re- 
fined level fields to the coarser levels and the interpolation of the boundary 
conditions for the refined levels from the coarser level information. Parti- 
cles at the boundary of the refined levels are regenerated according to the 
particle distribution at the coarser level. The details about inter-domain 
communications are provided in Sec. [5] 

• the method used to advance particles and fields on each domain is the 
Implicit Moment PIC method. 

The computational cycle involving collectively all the cloud of domains can 
then be represented as in Fig. |4j which depicts the modes and timing of the 
inter-level communications superimposed to the Implicit Moment method cycle 
of Fig. [I 

The computation starts on each domain with the moments equations being 
computed from the particle moments previously calculated. At this point, field 
information from the coarser levels is needed to provide the boundary conditions 
for the electric and magnetic fields on the refined grids ("Boundary Conditions 
E, B" step in the picture). After computation of the fields at time n + 1, the 
refined grid fields are projected to the coarser levels ("Refined Fields E, B"), thus 
updating their values with the increased level of accuracy granted by the refined 
grids. Then the computation resumes on each domain independently, with the 
interpolation of the field data to the particles, followed by the movement of the 
particles in the newly computed fields. When the new particle data are available 
at the coarser levels, particles are repopulated at the refined levels according to 
the parent distributions ("Particle repopulation" ) , thus enabling the calculation 
of particle moments on all levels. Notice that the repopulation of particles at the 
refined levels must be done before the calculation of particle moments, since the 
already updated density must be used to enforce Gauss's law before advancing 
the fields at the next time step. 

5. Inter-domain interaction 

In this Section, the inter-domain communications represented by vertical 
lines in Fig. |4] are explained in details. Three communication steps are needed 
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Figure 4: Representation of one computational cycle advancing particles and 
fields on each domain of the cloud. Each domain is completing the same opera- 
tions as all the others, with three points along the pipeline of operations (marked 
by vertical arrows) where communications across the levels become necessary. 
The fields from the coarser grids are used to provide boundary conditions for the 
refined grid fields ("Boundary Conditions E, B") before the field solver starts at 
the refined grid levels. After the field solution is known on the refined grids, it 
is projected to the coarser domains ("Refined Fields E, B"). After the particle 
mover step and before calculating particle moments at the refined grid level, 
the refined particles at the boundaries are regenerated according to the coarser 
grids particle distributions ("Particle repopulation" ) . 



across the grids to perform the following operations: 



project the fields from the refined to the coarser levels (Sec. 5.1 1; 



interpolate the boundary conditions for the refined levels electric and mag- 



netic fields from the coarser levels fields (Sec. 5.2 1 



• generate the refined level boundary particles from the coarser levels par- 
ticle distribution. Particles can leave the domain in the refined levels, but 
similarly new particles, created using information from the coarser lev- 
els, should be created in the boundary cells with positions and velocities 
derived from the coarser level particle distribution. Notice that particle 
information is exchanged only from the coarser to the refined grids (no par- 
ticle coalescence is needed, since particles are simulated also at the coarser 
grid levels) and that the particle moments (densities, currents, pressures) 
in the coarser grids are calculated natively from the local particles and 
not projected from the refined grids. 
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5.1. Projection from refined to coarser levels 

Information is transferred from the refined to the coarser grids through the 
"upwards" projection of the electric field and the consequent re-calculation of 
the magnetic field. 

The three components of the electric field are projected from the refined to the 
coarser grids as in 

where the subscripts gi and gi+i refer to the coarser and refined grid respec- 
tively, the coarsest grid being labeled with the index 0, Ep denotes the projected 
electric field and Ejv the electric field calculated natively on the grid. Notice 
that Ep_gj on the coarser grid is obtained as the average of two components: the 
electric field Eat ^j calculated natively on the coarser grid and the projection to 
the level gi of the electric field calculated at the level gi+i- Such projection is 
done as if each node in the refined grids was a particle in the coarser grids, with 
a normalizing factor as denominator: Wg^ is the same interpolation function 
defined in Eq. [12] and used in Eq. [TO] to obtain the particle sources on the grid 
points starting from the positions and velocities of the particles. 



The projection method of Eq. 30 is a key difference between the MLMD ap- 
proach presented here and the more common AMR algorithms which discard 
the information from the coarser levels: here the fact that both the coarser and 
the refined levels are fully functional allows to combine the information from 
both instead of simply neglecting the coarser grid information. 
Notice that a similar approach has been adopted, albeit only in a limited por- 
tion of the coarse domain, in Ref. |37j . where a coarse grid simulated with 
Magneto-Hydro-Dynamic (MHD) equations is interlocked with a refined PIC 
domain, as part of a strategy to eliminate high-frequency noise injected from 
the refined to the coarser domain through the upwards projection of the fields. 
Indeed, it is similarly observed here that calculating Ep^^, as an average instead 
of simply substituting it is beneficial to the evolution of the MLMD system in 
cases where the growth of a strong longitudinal electric field is not expected 
(e.g.: the Maxwellian, Weibel and shock benchmarks described in Sec. [6]). In 
fact, when the projected electric field is calculated by substitution, the particle 
density rig^ and the native field Ejv,g, in the coarser grid area which overlaps the 
refined grid develop a significantly higher level of noise when compared to the 
neighboring cells, thus introducing the risk of using non optimal electric field 
values as boundary conditions for the refined grids. It is easy to understand 
why a noisier density rig, degrades E^"*"^^ in the area of overlap: the divergence 
cleaning step |33j uses a degraded rig^ to " correct" Ep , which is later used to 
compute the Right Hand Side of Eq 



26 



and calculate the native field E^. at 

time n - 

It is less straightforward to understand why, when projecting the electric field 
by substitution, the particle density rt^, is noisier in the overlap area, since 
particles are moved in the projected and not in the native fields. A possible 
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explanation is that, when the simulated problem does not cause the growth of a 
strong longitudinal field, the longitudinal field on the coarser and refined levels 
reduces to noise fluctuations whose sign and amplitude may be different for the 
coarser and refined levels at a corresponding coordinate. This may introduce 
spurious effects in the coarser level particle motion, since particles are moved 
with the native coarser grid field in the area of no overlap and with the pro- 
jected refined grid field in the overlap area: consistency issues may arise at the 
boundary area, when particles are subject to the abrupt change from the native 
to the projected field and vice versa, and be amplified during the following steps 
of the calculation. 

When the electric field projection is done by averaging, instead, the particle 
density and the native electric field in the overlap area exhibit a level of noise 
comparable to the one in the neighboring areas. This is due to the fact that half 
of the native field information is retained when calculating the projected field 
and thus an effective coupling between the refined and coarser grid is reached, 
allowing particles to experience a smoother field transition when moving across 
the overlap area. 

The magnetic field in the overlap area of the coarser grids is calculated accord- 



ing to the discretized Maxwell-Faraday equation (Eq. 21 1) from the projected 
electric field: 

BJ+;=B?,,^,-AiVxE^+^ (31) 

The calculation of the magnetic field through Eq. [3l] is preferred to the most 
immediate alternative, an average between the coarser and refined grid fields 
similar to Eq. [30j since the Maxwell- Faraday equation at grid level is considered 
paramount to preserve. 

As regards the timing of the projection operations, it is important to remind 
(see also Fig. |4]) that they are carried out before moving particles in the coarser 
grid in order to increase the consistency between projected fields and particles 
which, as it was argumented before, is fundamental in the MLMD system. The 
base equations for particle motion on the coarser grids, originally expressed as 



in Eq. 17 thus become 



•p -p 



TO., 



(E';+;(x;'+V2)+,n+l/2^B?,,^(^„+l/ 



(32) 



where the projected fields Ep^^^ and Bp^^ are used instead of the native 
ones in a Predictor Corrector approach. The finest grid particles are moved, of 
course, in the native fields. 

5.2. Interpolation from coarser to refined levels 

In order to drive the evolution of the refined levels, their boundary points 
are obtained by interpolation from the fields of the level above. Again the grid 
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points of the refined levels are treated by the coarser levels as if they were 
particles and the interpolation is therefore: 

V'/.gi + i = "^^^91 ^91 (Xff, - ^91 + 1 ) (33) 

9l 

where V'/.gi+i denotes the interpolated point on the refined grid and "0 is a generic 
notation for both the electric and magnetic field. Notice that the boundary con- 
ditions for the refined grid are calculated from the native, not projected, coarser 
grid fields '^n.qi and that the points needing interpolation are the boundary 
nodes, e.g., the ghost cells nodes. This minimal information exchange for fields, 
limited to the ghost cells nodes, is sufficient to drive the refined grid evolution 
also in challenging cases as the shock benchmark presented in Sec. |6.4[ 
Since the boundary conditions for the electric field are needed at the refined 
grid level before calculating E^"*"^^^^ , a bottleneck for a future, planned parallel 
evolution of the presented algorithm is here foreseeable: the refined grids have 
to wait for the calculation of the coarser grid 'E^g^ before starting their own 
calculation. This bottleneck can be overcome by using as boundary conditions 
for the refined grid not the final result of the iterative calculation of Eq. [26] on 
the coarser grid, but an intermediate one. The assessment of the viability of 
this solution is left as future research. 

5.3. Particle repopulation in the refined levels 

In the presented MLMD approach, particles information has to be exchanged 
between the coarser and the refined grids only in the "downwards" direction: 
since the coarser levels are fully simulated, particles exiting the refined domains 
can be simply lost, without the need of using them to recreate particles on the 
coarser grids. However, refined grid particles have to be recreated at the bound- 
aries of the grids consistently with the coarser grids particle distribution. 
Two possibilities are open for particle repopulation: particles can either be 
recreated according to a fixed distribution function, e.g. a Maxwellian distri- 
bution function, with parameters (drift velocity, thermal spread) derived from 
the coarser grid, as in Ref. |37] and Ref. [35] (notice however that in both these 
cases the coarser level is not simulated with a particle method and therefore full 
particle information is not available), or by taking full advantage of the fact that 
the velocities and positions of coarser grid particles are known, as in Ref. |7]. 
In taking this decision, two factors have been considered: as a first attempt, 
we decided to retain the possibility to simulate non Maxwellian distribution 
functions at the grid interface, option which is lost if the first approach is 
adopted [35] • Secondly, we considered that particle repopulation is a task which 
requires particular attention since consistency in particle motion between the 
refined and coarser grids is not enforced explicitly in the MLMD algorithm pre- 
sented here, but comes as a consequence of the projection and interpolation 
operations described in Sec. |5.1| and Sec. |5.2[ 

For this reason, a rather conservative approach has been opted for in regener- 
ating particles in order to guarantee the consistency between the particle pop- 
ulation on the refined and coarser grids at least at the grid boundaries (further 
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considerations about particle motion consistency across the grids will follow in 



Sec. 6.3) 



First of all, particles are regenerated not only in the ghost cells, but also in a 
small number of neighboring cells already inside the active part of the refined 
grid. The total area where particles are regenerated is called here " Particle Re- 
population Area" (PRA). This operation, which may seem redundant, proves 
instead to be of fundamental importance in cases when boundary particles are 
subject to the formation of structures in phase space, as in the Weibel and two 
stream instability benchmarks shown in Sec. [6j The particles populating the 
PRA from the previous time step are deleted and new particles are generated 
by reproducing the parent particle distribution in the corresponding area of the 
parent grid. Each particle in the parent grid sitting in an interval corresponding 
to the PRA ± dxg^, where dxg^ is the coarser grid spacing, is split, on the re- 
fined grid, into RF particles {RF being the Refinement Factor, that is the ratio 
between the coarser and the refined grid dimensions) following the algorithm 
described in Ref . [2D] , used in Ref . [H] and summarized here in Eq. 34p6 for the 



ID case (the extension to higher dimensions is straightforward, just requiring 
to adapt Eq. [36] in order for the combined shape function of the refined grid 
particles to reproduce the coarser grid particle shape function): 

Ct,\=CV^^ (34) 



Psi + i Pgi 



(35) 



X 



■ dxg^^-^ + ^ ^0,1+1, for i = : RF — 1 (36) 



Eq. [34] guarantees that the total charge is respected, Eq. 



35] that the velocity 
distribution is reproduced without distortions and Eq. 36 where xo.j+i is the 
coordinate referred to the coarser grid at which the active part of the refined 
grid starts, that the combined shape function of all the RF refined grid particles 
corresponds to the shape function of the coarser grid parent particle. 
Notice that, as visible in Fig. [4] and remarked by the temporal index n -I- 1 in 
Eq. |34p6[ the particles in the refined grid are generated from the parent grid 
distribution after the particle mover operations have been completed on both 
grids. The aim is to calculate particle moments consistent with the repopulated 
particle distribution on the refined grids. 

As regards the number of cells to repopulate, that is, the size of the PRA 
area, the rule of thumb followed in Sec. [6] is to dimension the repopulation 
area according to the distance that the particle with the maximum reasonably 
expected velocity in the system can cover in a time step. Future implementations 
will make the choice of the number of PRA cells automatic. 
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6. Benchmarks 



The simulation of a MLMD PIC system entails a variety of challenges, which 
span from checking fields and particle continuity at the grid boundaries to the 
more delicate tasks of assessing the impact of the exchange of information be- 
tween the levels, verifying the results of the simulations against theoretical ex- 
pectations and understanding which degree of consistency in the evolution of 
the different grid levels should be expected. 

These issues are addressed here through a series of simulation challenges in a 
ID environment, each devoted to test a specific aspect of the system: 

• how the system reacts to the exchange of information between the levels. 
A Maxwellian plasma is simulated in a MLMD environment to examine 
how the two coupled grids react to the exchange of information in absence 
of notable plasma activity; 

• how an instability is simulated across two levels and how not moving 
particle structures are created and preserved across the grid boundaries. 
A Weibel instability is excited in the coarser and refined grids and its 
evolution is followed. In particular, the growth rate of the instability and 
the continuity of fields and particles across the grids are checked, with 
particular attention to the development of electron structures in phase 
space; 

• how moving particle structures react to crossing the grid boundaries. A 
two stream instability is simulated and the continuity of electron holes 
moving across the coarser-refined grid boundary is checked. Moreover, 
the comparison of the electric field and electron holes evolution of the two 
grids will stir considerations about particle motion consistency in the two 
grids; 

• how the refined grid reacts to strong driving from the coarser grid. A 
shock is excited at the coarser grid boundaries and launched towards the 
center of the coarser grid, where the refined grid sits. The propagation of 
the shock wave across the domains is studied. 



6.1. Testing the system reaction to the exchange of information between the 
refined and coarser grids: Maxwellian plasma simulation 

To test how grid coupling influences the simulation evolution in absence 
of notable plasma activity, a Maxwellian plasma is simulated across the two 
domains. Electrons and ions with a realistic mass ratio of mi/me = 1836, a 
thermal velocity of Vth/c = 0.2 in all directions (c is the speed of light, used as 
normalization factor for the velocities) and no drift are immersed in a domain 
of size Lx,gio/de — 84, with de the electron skin depth. A refined grid with 
Refinement Factor RF=4 overlaps the coarser grid at 21 < x/d^ < 42. The 
information exchange between the two grid is done accordingly to Sec.[5j with a 
PRA (see Sec. 5.3 1 extending for one cell into the active part of the refined grid. 
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Both grids have 280 cells, the time step is ujpedt = 0.15, with ujpe the electron 
plasma frequency, and the simulation is carried on for 3000 cycles. A reference 
case with two levels simulated independently, both with periodic conditions for 
fields and particles, is also shown to help assessing the impact of the MLMD 
algorithm on this basic simulation. 

Fig. [5] depicts the total field energy normalized to the initial energy (hence the 
low value: the energy of the system is mostly stored in particles) for the MLMD 
(panel a) and the reference (panel b) case respectively, with the blue line referred 
to the coarser and the red line to the refined grid in both panels. 
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Figure 5: Effect of the Multi Level Multi Domain system on the total field 
energy in the simulation of a Maxwellian plasma. Total field energy normalized 
to the total domain (fields plus particles) energy at Upet = for the (a) Multi 
Level Multi Domain and (b) reference case respectively. The blue line is relative 
to the coarser grid, the red line to the refined grid. 



Notice that in the reference case, where no interlevel communications are 
present (panel b), the field energy stabilizes around different values for the 
coarser and refined grid, with the coarser grid exhibiting higher energy values. 
Instead, in the MLMD case (panel a) the two grids tend to converge to the 
coarser grid energy level after a transient time which corresponds to the time 
needed for the repopulated particles to overcome in number the native ones 
which exit the grid by thermal motion. The fact that the refined level energy 
raising to the coarser grid values is connected to particle repopulation is con- 
firmed also by check simulations with MLMD conditions on the refined grid 
fields but periodic conditions on the refined grid particles, which exhibit energy 
levels comparable to the reference case. 

This evolution of the field energy in the MLMD system is considered a symp- 
tom of the efficient coupling between the refined and coarser grid, which is 
mostly achieved through particle repopulation. A spectral analysis of Ex^N.gn 
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and Ex^p^gig in the MLMD system, additionally, shows that this effect is not 
detrimental to the physical significance of the simulation. 

Fig. |6] panel a, shows the Fast Fourier Transform of E^^N.gn in the MLMD 
system, while Fig. [6j panel b, depicts the difference of it with the corresponding 
spectrum for the reference case. 
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Figure 6: Effect of the Multi Level Multi Domain system on the spectrum of 
the longitudinal electric field in the simulation of a Maxwellian plasma, (a) Fast 
Fourier Transform of Ex.N,gii in the Multi Level Multi Domain simulation and 
(b) difference with the corresponding spectrum in the reference case. 

Notice that Fig. [6^ exhibits all the features expected from an implicit PIC 
simulation of a Maxwellian plasma: the Langmuir wave (marked as (1)) and 
the dispersion relation for the light wave (marked as (2)) as reproduced by an 
implicit PIC code|26| ^marked as (3)) are both visible and not modified in the 
MLMD case when compared to the reference one, as shown in Fig. [6}d. The 
extra energy in the refined grid of the MLMD case is stored at low k/deS in an 
area unaffected by relevant features and thus proves not to damage the physical 
relevance of the simulation. 
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A similar analysis carried out for the coarser level shows that no relevant differ- 
ence is observable between the spectrum of Ex.p.g^g in the MLMD and reference 
case. 

6.2. Testing the development of instabilities and the continuity of stable particle 
structures across the grid boundary: Weibel instability simulation 
The Weibel instability is excited by an anisotropy in the thermal velocity of 
particles, i.e. by the presence of a preferential, warmer direction in the parti- 
cle distribution function. The resulting particle motion produces a current in 
the direction of higher temperature which in turn grows, through the Ampere's 
equation, a perpendicular magnetic field. The instability eventually saturates 
due to magnetic trapping processes which lead to the formation of characteristic 
structures in the phase spaces of electrons, with trapping points being formed at 
the zeros of the magnetic field and a consequent zig-zag distribution emerging 
in the Vy/c vs x/d^ electron phase space [521 EHl SI]- In a ID system simulated 
in the x direction with y the direction of higher particle temperature, a current 
Jy is formed in the warmer y direction and a magnetic field Bz is consequently 
grown. 

It is investigated here how the magnetic field is originated in the MLMD system 
and the continuity of the particle structures across the grid boundaries. 
The Weibel instability is excited by setting the thermal velocities for electrons 
to Vth,y/c — 0.2 and Vth,x/c — Vth,z/c = 0.05 in the y and x and z direction 
respectively. The ions, with a mass ratio of mi /me = 1836, are loaded with the 
same temperature and temperature anisotropy as the electrons. The time step is 
LUpedt = 0.1, the coarser grid has length L^.g^g/de — 18.47 and both the coarser 
and the refined grids have 280 cells. The refined grid overlaps the coarser grid 
at 6.13 < x/de < 10.75, the Refinement Factor is RF — 4 and three PRA cells 
are used in the active domain of the refined grid. Remember that the number of 
PRA cells to be used is decided a priori according to the spatial and temporal 
resolution of the refined grid and the expected maximum velocity of particles. 
More sophisticated algorithms will be devised for future implementations. 
From the linear theory, the dominant wavenumber in such a configuration is 
kde = 1.36 with growth rate j/ujpe = 0.123. In the presented simulation geom- 
etry, the dominant mode n, with L^.g^, — n'^, is n = 4. 

Fig. [t] shows the current Jy (panel a) and the magnetic field (panel b) in 
the aforementioned MLMD case, with the refined grid fields superimposed to 
the coarser grid ones. Notice that both Jy and B^ develop as expected [iU] . 
with a calculated growth rate for kde = 1.36 equal to j/ujpe = 0.123, the the- 
oretical value, and that no discontinuity is noticeable at the boundaries of the 
refined grid. It is therefore possible to state that the correct development of the 
instability is preserved by the MLMD algorithm. 

More interesting is however the inspection of the phase space plots, with the 
aim of checking if the electron trapping structures are correctly developed and 
how they are dealt with across the grid boundaries. 

Fig. [s] and |9] show, in panel al and all respectively, the phase space v^/c vs 
x/de and Vy/c vs x/c?e for the electrons of the coarser grid alone, while in panel 
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Figure 7: Relevant physical quantities in the Multi Level Multi Domain simu- 
lation of the Weibel instability, (a) Current Jj/,g,o with superimposed Jy,g„ and 
(b) magnetic field -B^.p,^,,, with superimposed Bz,N,gn at 6.13 < x/de < 10.75. 



bl and bll the refined grid phase spaces are superimposed to the coarser grid 
ones to check for particle structure continuity. Fig. [8] refers to time uipet = 32, 
when particle trapping starts becoming evident, while Fig. [9] is taken after the 
saturation of the mode with n = 4, at ujpet = 48 (saturation occurs at ujpet = 41). 
Note that the particle structures develop as expected from previous studies [ICTl 
I41| and as in the reference case, not shown here, with electron trapping points 
being formed at the zeros of in the simulated x direction and a zig zag 
distribution emerging in Vy/c vs x/de in accordance with the current structure. 
Notice that the particle structures are coherent between the refined and the 
coarser level and continuous among the refined and coarser grid boundaries. 
The Vy/c vs x/de phase space plot in Fig. [9] panel bll reveals a particularly 
impressive grid coupling: even the smallest scale electron features are correctly 
captured across the grid boundaries. 
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6.3. Testing the continuity of moving particle structures across the grid bound- 
ary: two stream instability simulation 
While in the case of the Weibel instabiUty the continuity across the refined 
and coarser grid boundaries of the electron structures in phase space could be 
argued to be a consequence of the fact that the trapping positions for particles 
are located at fixed points, similar arguments cannot be made in the case of 
a MLMD simulation of a two stream instability ^42^. In this case in fact the 
evolution and merging of the generated electron holes cannot be determined a 
priori and the structures move through the domain, since the instability has a 
real frequency. 

In the case here presented, the two stream instability is excited through an 
electron stream velocity of Vx/c — ±0.2 in the simulated direction, while the 
thermal velocities for ions, with realistic mass ratio, and electrons are Vth/c = 
0.04. The coarser grid has length L^^g^g/de = 11.94, the time step is ojpedt — 0.1 
and 2000 computational cycles are executed. The refined grid has Refinement 
Factor RF = 4, extends at 3.96 < x/d^ < 6.95 and nine PRA cells are used 
(such an high number of PRA cells is due to the high velocities expected in the 
system after the development of the instability). Fig. 10 depicts the electric field 
Ex,p,gio the coarser grid in panel a, while in panel b the refined grid Ex.N,gii 
is superimposed at 3.96 < x/de < 6.95. In both cases, the contour of the coarser 
field density Ug^g is shown on top of the electric field. This picture, contrarily to 
before, does not depict the entire coarse grid domain and simulation duration, 
but just a fraction of both: 3.96 — 0.5 < x/de < 6.95 + 0.5, with a ±0.5/de 
distance from the finer grid boundary, and LUpet < 80, when the electron holes 
merging activity is more intense. The aim is to show that while, as usual, field 
continuity is granted at the boundary between the coarser and refined grid, a 
certain mismatch in is visible inside the overlap area. Focus, for example, 
on the area highlighted in red in both panels: notice that the traces of negative 
electric field, even if roughly similar, exhibit small differences. 

This mismatch is due to the fact that the electron holes that produce the 
traces in evolve slightly differently in the two grids, since the coarser grid 
forcing over the refined grid is limited to the boundaries for fields and to the 
few PRA cells for particles. Therefore, some differences in the particle evolu- 
tion and, consequently, in the traces in E^ in the overlap area far away from the 
boundaries have to be expected. The above mentioned mismatch between the 
electric field evolution at different levels is removed if the electric field is pro- 
jected by substitution instead of averaging. However, as already argumented in 
Sec. |5.1[ projecting the electric field by substitution proves detrimental in cases 
where the major simulated dynamics are not electrostatic. For this reason, for 
sake of generality of implementation, the projection of the electric field is done 
by average as described in Sec. |5.H with the caveat that particle moments and 
fields are consistent within the grids (i.e., Epg^ is consistent with the density 
rigj on the same level) and not across the grids (i.e., £^Ar,g,^.i is not necessarily 
consistent with the density rig, in the overlap area of the coarser grid) . This is 



well visible in Fig. 10 where the coarser grid density Ug^^ is superimposed to 



the coarser grid field in panel a and to the refined grid field in panel b. It is well 
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known that the peak density for electron holes is registered at the center of the 
electric field bipolar trace[lS]: see, focusing for example on the area enlightened 
by the red rectangle, that Ugig is consistent with the coarser field Ex^p^gif^, not 
with the refined grid one, jv.gi+i, even if the differences are really minimal to 
catch, ng,^, of course, is consistent with the refined grid fields. 
Let us check now what this means in term of phase space plots. 



Fig. 11 and Fig. 12 depict the phase space plots Vx/cvs x/de (marked as I) and 
Vy/c vs x/de (marked as II) for the above mentioned simulation at two different 
times, LUpet = 62 and ujpet = 92. In both cases, panel a refers to the phase 
spaces of the coarser grid, while in panel b the phase spaces for the refined grid 
are superimposed to the coarser grid ones, shown in transparency. 

Fig. 11 is taken at ujpet — 62, a time at which the differences between Ex^p^io 



and Ex^N,n are rather pronounced, how it is possible to check by following the 
red vertical line in Fig. [lOl panel a and b. Notice indeed that, comparing the 



Vx/cvs x/de phase space plots in Fig. 11 al and bl, it is possible to see the cause 



of the field mismatches at 5 < x/de < 6.1: the inner part of the central electron 
hole in the overlapping area is moving rightwards with an higher velocity in 
the coarser grid than in the refined grid, thus producing a shift in the electric 
field trace. Notice, however, that particles are continuous at the boundaries, as 
expected. 



Fig. 12 instead, is taken at a later time, ujpet = 92, when no mismatch in 
the electric field is visible between the grids. In this case, it is possible to fully 
appreciate the advantage that the MLMD simulation offers in terms of resolution 
over a simulation with the coarse resolution alone. Notice the increased level of 
details provided by the refined grid phase spaces and how the smallest streams 
at high, negative v^/c nicely cross the interface, both from the coarser to the 
refined grid and vice versa, in Fig. [T2j3l. 



6.4- Testing the reaction of the refined grid to strong driving from the coarser 
grid: shock simulation 
The test problems presented in Sec. |6.2| and |6.3| focused on the development 
of instabilities excited simultaneously on the refined and coarser grids: the issues 
under investigation were the consistency of the evolution on the two grids and 
the continuity of particle structures at the grid boundaries. 
In this case, instead, a shock problem is used to test the effectiveness of the 
coarser level driving over the refined grid: a shock wave is created at the coarser 
grid boundaries and then launched towards the center of the domain and thus 
the refined grid. 

The plasma inside the domain is initially in a Harris state [44 , with half width 
of the current sheet Lu/di = 0.5 {di = c/ujpi is the ion inertial length and Wpj 
the ion plasma frequency), a temperature ratio of Ti/Te = 5 and a mass ratio of 
nii/me — 25. The mass and temperature ratio are the same as in the Newton 
challenge [45 , while the half width of the current sheet is thinner in order to 
have a significative particle density only at the center of the coarse domain: the 
aim is to have the wave excited from the boundary compression (an electric 
field of arbitrary value Ey — 0.015 is imposed at the boundaries of the coarser 
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grid) propagating across the domain with hght speed. For the same reason, a 
particle background, which would decelerate the wave propagation velocity to 
V = E X B/B^, is not included. 

The MLMD simulation is performed with a coarser grid length of L^ g^^Jdi — 20, 
128 cells for both the coarser and the refined grids, a time step ujpidt = 0.1, a 
Refinement Factor RF = 4 and no PRA cells in the refined grid. The refined 
domain is situated at the center of the coarser grid, at 7.5 < x/di < 12.5. 
The MLMD evolution is compared with check runs performed with the same 
grid length and with the coarser and refined level resolution (128 and 512 cells, 
same grid length) respectively, all the other parameters unchanged. Notice that 
references runs were not shown (however performed) in the previous cases since 
the evolution of the systems presented previously is well known and documented. 



Fig. 13 depicts the evolution of Ey in the three cases, the coarse (panel a) 
and fine (panel b) reference cases and the MLMD system (panel c), where 
the refined grid field Ey^j^^g^^ is superimposed to the coarser grid field Ey^p^g^ 
between 7.5 < x/di < 12.5. Notice that in all three cases the wave propagates 
undisturbed towards the center of the coarse domain and is then reflected. No 
spurious reflections or front distortion are observable in the MLMD case. 



Fig. 14 shows the Jy current evolution again in the three cases, not refined 
(panel ajand refined (panel b) reference case and MLMD system (panel c). Only 
the area at 7 < x/di < 13 is shown, since no significative activity is detectable 
further away from the center of the domain. In all three cases it is possible to 
observe how the arrival of the shock modifies the current profile starting from 
ujpit = 10, with an initial thicking of the current profile and a consequent peak. 
Notice how the MLMD simulation allows to reach the same resolution level for 
the current as in the refined case in the high density area of the simulation, 
while removing the need of resolving to an unnecessary level the less relevant 
areas of the simulation. 



7. Conclusions 

In this paper, a novel solution of the problem of simulating larger systems 
for longer times and with kinetic scale resolutions is proposed. 
The novelty of the approach presented is twofold. First, a Multi Level Multi 
Domain solution is introduced: the simulated levels are fully functional and self- 
similar and, even when a more refined grid exists for an area, the corresponding 
coarser area is still simulated completely, with both fields and particles (see 



Sec. 3.1). Second, the framework in which the MLMD technique described 



here is applied is the Implicit Moment PIC method (see Sec. 3.2 and references 
therein for a review of the method) . 

The choice of adopting a MLMD instead of an AMR approach is one of the 



possible solutions (see Sec. 2.1) to the problem of adapting the particle shape 
size to the grid size. Moreover, having fully functional self-similar levels is a 
great advantage from the point of view of the ease of programming (with an 
object-oriented approach, each domain is treated as an instance of the same ob- 
ject) and remarkably simplifies the daunting task of managing the boundaries 
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between the coarser and refined grids. 

Implicit PIC methods already grant a computational advantage over explicit 
methods since they notably relax the stability-dictated limits on grid resolu- 
tions: they are therefore optimal candidates for AMR and Multi Level evolu- 
tions, given the possibility of spanning a notably increased range of time and 
space steps. 

The algorithm proposed in the paper prescribes an exchange of information be- 
tween the coarser and refined grids in three steps (see Sec. [s]): projection of 
the refined grid fields from the refined to the coarser grid, interpolation of the 
refined grid boundaries from the coarser grid fields and particle rcpopulation at 
the refined grid boundaries starting from the particle distribution in the coarser 
grid. 

The key point in exchanging field information between the grids (i.e., both 
for projection and interpolation operations) is to consider the grid points in 
the refined grids as particles in the coarser grid and consequently use for both 
projection and interpolation the same interpolation functions employed for cal- 
culating the momenta from particle positions and velocities. 
The projection strategies presented have the primary aim to favor the maxi- 
mum coupling between the levels while preserving fields and particles consis- 
tency within the grids. For this reason, the electric fields are projected from 
the refined to the coarser grids as an average between the native, coarser fields 
and the refined grid information (Eq. 30 1 and the magnetic fields are recalcu- 
lated from the projected electric field through the Maxwell-Faraday's equation 



(Eq. 31). Eq. 30 makes the transition that coarser level particles experience 



when crossing into the overlap area smoother, while Eq. |3T] is adopted in or- 
der to safeguard the validity within the grid of the Maxwell-Faraday equation 
iEq.§. 

The strategy adopted in particle repopulation, then, based on the particle re- 
zoning algorithm of Ref. |20) . is chosen to favor the grid coupling by exactly 
reproducing the coarser grid particle distributions at the boundaries of the re- 
fined grids. 

The proposed algorithm is tested in a ID setting against a series of challenges. 



In Sec. 6.1 the reaction of the MLMD system to the exchange of information 
between different levels in absence of plasma activity is tested through the sim- 
ulation of a Maxwellian plasma. Sec. |6.2| and |6.3| focus on the development of 
instabilities and the formation of particle structures either fixed in position but 
grown across the refined-coarser grid boundaries (MLMD simulation of a Weibel 
instability) or moving across the grid boundaries (MLMD simulation of a two 
stream instability). Finally, in Sec. 6.4 a shock is excited in the coarser grid 
and launched across the finer grid to test the effectiveness of the driving at the 
boundary. 

The results are encouraging: the proposed MLMD algorithm proves to grant 
an excellent coupling between refined and coarser levels (Sec. 6.1) and not to 
impair the development of instabilities and particle structures across the grids 



(Sec. 6.2 and 6.3). Moreover, the shock test shows that even an extreme case. 



like a wave propagating at light speed across the coarser-refined grid boundaries. 
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is correctly handled by the system and, even more relevantly, that the aim of 
the MLMD system is fulfilled. In cases when only a small section of the simula- 
tion requires enhanced resolution, while the rest may allow the use of a coarser 
grid, the algorithm proposed proves capable of correctly driving the refined grid 
towards the evolution expected in a simulation fully performed with the refined 
level resolutions, without wasting computational resources in the areas which 
requires less resolution. 

Future work will focus on the application of the proposed algorithm to the op- 
erative implicit PIC code Parsek2D [29] and on the quantitative assessments of 
the performance increase granted by the MLMD method when used in a parallel 
implementation. 



Acknowledgments 

The present work is supported by the Exascience Intel Lab Europe, by the 
Onderzoekfonds KU Leuven (Research Fund KU Leuven) and by the Euro- 
pean Commission's Seventh Framework Programme (FP7/2007-2013) under the 
grant agreement no. 263340 (SWIFF project, www.swiff.eu) and No. 218816 
(SOTERIA project, www.soteria-space.eu) 



References 



W. Gonzalez, B. Tsurutani, A. Cla de Gonzalez, Interplanetary ori- 
gin of geomagnetic storms. Space Sci. Rev. 88 (1999) 529-562, doi: 
10.1023/A:1005160129098. 

A. S. S. Lipatov, The Hybrid Multiscale Simulation Technology, Springer, 
Berlin, 2002. 

D. Krauss-Varban, N. Omidi, Large-scale hybrid simulations of the mag- 
netotail during reconnection, Geophys. Res. Lett. 22 (1995) 3271-32-74. 

C. Birdsall, A. Langdon, Plasma Physics Via Computer Simulation, Taylor 
& Francis, London, 2004. 

R. Hockney, J. Eastwood, Computer simulation using particles, Taylor & 
Francis, 1988. 

Giovanni, Lapenta, Particle simulations of space weather, J. Comput. Phys. 
231 (3) (2012) 795 - 821. |doi : 10 . 1016/ j . j cp . 2011 . 03 . 035) 

K. Fujimoto, R. D. Sydora, Electromagnetic particle-in-cell simulations 
on magnetic reconnection with adaptive mesh refinement, Comput. Phys. 
Commun. 178 (2008) 915-923. doi : 10 . 1016/j . cpc . 2008 . 02 . 010 

J. U. Brackbill, An adaptive grid with directional control, J. Comput. Phys. 
108 (1993) 38. 



32 



doi : 



[9] J. Brackbill, An adaptive grid with directional control for toroidal plasma 
simulation, in: Proc. 14th International Conference on the Numerical Sim- 
ulation of Plasmas, American Physical Society, Annapolis MD, 1991, 1991. 

[10] G. Lapenta, Democritus: An adaptive particle in cell (pic) code for object- 
plasma interactions, J. Comput. Phys. 230 (12) (2011) 4679 - 4695. 
ri0.1016/j -jcp. 2011 .02 .0411 

[11] M. J. Berger, J. Oliger, The AMR technique, J. Comput. Phys. 53 (1984) 
484-512. 

[12] J. Vay, P. Colella, A. Friedman, D. P. Grote, P. McCorquodale, D. B. 
Serafini, Implementations of mesh refinement schemes for Particle-In-Cell 
plasma simulations, Comput. Phys. Commun. 164 (2004) 297-305. 
il0.1016 /j . cpc . 2004 . 06 . 075 

[13] G. B. Jacobs, J. S. Hesthaven, High-order nodal discontinuous Galerkin 
particle-in-cell method on unstructured grids, J. Comput. Phys. 214 (2006) 



doi : 



96-121. doi : 10 ■ 1016/ j . j cp . 2005 . 09 . 008 



[14] G. Lapenta, J. U. Brackbill, P. Ricci, Kinetic approach to microscopic- 
macroscopic coupling in space and laboratory plasmas, Phys. Plasmas 13 
(2006) 055904. 

[15] Y. N. Grigoryev, V. A. Vshivkov, M. P. Fedoruk, Numerical particle-in-cell 
methods: theory and applications, VSP BV, AH Zeist, 2005. 

[16] J. Vay, P. Colella, P. McCorquodale, B. van Straalen, A. Friedman, D. P. 
Grote, Mesh refinement for particle-in-cell plasma simulations: Applica- 
tions to and benefits for heavy ion fusion. Laser. Part. Beams 20 (2002) 
569-575. Idoi : 10 . 1017/S0263034602204139, 

[17] P. Colella, P. C. Norgaard, Controlling self- force errors at refinement 



boundaries for AM R-PIC, J. Comput. Phys. 229 (2010) 947-957. [d^ 
|10 ■ 1016/ j ■ j cp ■ 2009 ■ 07 . 004[ 

[18] C. D. Boor, A practical guide to splines, Springer, 1978. 

[19] G. Lapenta, Variational grid adaptation based on the minimization of local 
truncation error: Time independent problems, J. Comput. Phys. 193 (2004) 
167-194. 

[20] G. Lapenta, Particle rezoning for multidimensional kinetic particle-in-cell 



simulations, J. Comput. Phys. 181 (2002) 317-337. |doi : 10 . 1006/jcph. 
l2002T7l26t 

[21] K. Fujimoto, S. Machida, Electromagnetic full particle code with adaptive 
mesh refinement technique: Application to the current sheet evolution, J. 



Comput. Phys. 214 (2) (2006) 550 - 566. |doi : 10 . 1016/ j . j cp . 2005 . 10 . 
10031 



33 



[22] R. Mason, Implicit moment particle simulation of plasmas, J. Comput. 
Phys. 41 (1981) 233. 

[23] J. Denavit, Time-filtering particle stimulations with uj^peAt ^ 1, J. Com- 
put. Phys. 42 (1981) 337. 

[24] J. Brackbill, D. Forslund, Simulation of low frequency, electromagnetic 
phenomena in plasmas, J. Comput. Phys. 46 (1982) 271. 

[25] A. Langdon, B. Cohen, A. Friedman, Direct implicit large time-step particle 
simulation of plasmas, J. Comput. Phys. 51 (1983) 107-138. 

[26] J. U. Brackbill, D. W. Forslund, An implicit method for electromagnetic 
plasma simulation in two dimensions, J. Comput. Phys. 46 (2) (1982) 271 
- 308. doi:D0I:10.1016/0021-9991(82)90016-X, 

[27] H. X. Vu, J. U. Brackbill, CELESTEID: An imphcit, fully-kinetic model for 
low- frequency, electromagnetic plasma simulation, Comput. Phys. Comm. 
69 (1992) 253. 

[28] S. Markidis, E. Camporeale, D. Burgess, Rizwan-Uddin, G. Lapenta, 
Parsek2D: An Implicit Parallel Particle-in-Cell Code, in: N. V. Pogorelov, 
E. Audit, P. Colella, & G. P. Zank (Ed.), Astronomical Society of the 
Pacific Conference Series, Vol. 406 of Astronomical Society of the Pacific 
Conference Series, 2009, pp. 237 — h. 

[29] S. Markidis, G. Lapenta, Rizwan-uddin, Multi-scale simulations of plasma 
with iPIC3D, Math. Comput. Simulat. 80 (7) (2010) 1509-1519. doi:DOI:j 
|lO. 1016/j .matcom. 2009 .08 .0381 

[30] H. Vu, J. Brackbill, Accurate numerical solution of charged particle motion 
in a magnetic field, J. Comput. Phys. 116 (1995) 384. 

[31] K. Noguchi, C. Tronci, G. Zuccaro, G. Lapenta, Formulation of the rel- 
ativistic moment implicit particle-in-cell method, Phys. Plasmas 14 (4) 
(2007) 042308— |doi : 10.1063/1 .2721083. 

[32] C. T. Kelley, Iterative methods for linear and nonlinear equations, SIAM, 
Philadelphia, 1995. 

[33] P. Ricci, G. Lapenta, J. Brackbill, A simplified implicit maxwell solver, J. 
Comput. Phys. 183 (2002) 117-141. |doi : 10 . 1006/ j cph . 2002 . 7170 

[34] Y. Saad, Iteractive Methods for Sparse Linear Systems, SIAM, Philadel- 
phia, 2003. 

[35] J. U. Brackbill, B. I. Cohen (Eds.), Multiple Time Scales, Academic Press, 
Orlando, 1985, Ch. Simulation of Low-Frequency, Electromagnetic Phe- 
nomena in Plasmas, (J.U. Brackbill, D.W. Forshmd), p. 271. 



34 



[36] Keizo, Fujimoto, A new electromagnetic particle-in-cell model with adap- 
tive mesh refinement for high-performance parallel computation, J. Corn- 
put. Phys. 230 (23) (2011) 8508-8526. doi : 10 . 1016/j . j cp . 2011 . 08 . 002, 

[37] T. Sugiyama, K. Kusano, Multi-scale plasma simulation by the interlocking 
of magnetohydrodynamic model and particle-in-cell kinetic model, J. Com- 



put. Phys. 227 (2) (2007) 1340 - 1352. doi : 10 . 1016/j . j cp . 2007 . 09 . Oil 



M. Shay, J. Drake, B. Dorland, Equation free projective integration: A 
multiscale method applied to a plasma ion acoustic wave, J. Comput. Phys. 



226 (1) (2007) 571 - 585. doi : 10 . 1016/j . j cp . 2007 . 04 . 016 



doi : 



[39] E. S. Weibel, Spontaneously growing transverse waves in a plasma due 
to an anisotropic velocity distribution, Phys. Rev. Lett. 2 (1959) 83-84. 
Fdoi : 10 ■ 1 103/PhysRevLet t . 2 . 83, 

[40] A. Stockem, M. E. Dieckmann, R. Schlickeiser, PIC simulations of the 
thermal anisotropy-driven weibel instability: field growth and phase space 
evolution upon saturation, Plasma Phys. Contr. F. 51 (2009) 075014. 
[10 ■ 1088/0741-3335/51/7/075014, 

[41] M. E. Innocenti, M. Lazar, S. Markidis, G. Lapenta, S. Poedts, Electron 
streams formation and secondary two stream instability onset in the post- 
saturation regime of the classical weibel instability, Phys. Plasmas 18 (2011) 
052104. doi : 10 . 1063/1 . 358708"9l 

[42] F.F.Chen, Introduction to Plasma Physics and Controlled Fusion, Plenum 
Press, 1984. 

[43] B. Eliasson, P. Shukla, The dynamics of electron and ion holes in a coUi- 
sionless plasma, Nonlinear Proc. Geoph. 12 (2005) 269. 

[44] E. Harris, On a plasma sheath separating regions of oppositely directed 
magnetic field, II Nuovo Cimento (1955-1965) 23 (1962) 115-121. 

[45] J. Birn, K. Galsgaard, M. Hesse, M. Hoshino, J. Huba, G. Lapenta, P. L. 
Pritchett, K. Schindler, L. Yin, J. Buchner, T. Neukirch, E. R. Priest, 
Forced magnetic reconnection, Geophys. Res. Lett. 32 (2005) L06105. 



35 



Grid 0, (0 t= 32 




2 4 6 



(bl) 



Superimposed, m t= 32 



0.2 
0.1 


-0.1 
-0.2 



mm 



2 4 6 



x/d 



10 12 14 16 18 



(bll) 



Superimposed, (0^^!= 32 



0.5 



-0.5 




2 4 6 



x/d 



10 12 14 16 1£ 



I 
I 



Figure 8: Phase spaces before saturation for the Multi Level Muhi Domain 
simulation of the Weibel instability, (a) Phase space Vx/c vs x/de (al) and Vy/c 
vs x/de (all) for the coarser grid and (b) phase space v^/c vs x/de (bl) and 
Vy/c vs x/de (bll) for the refined grid superimposed to the coarser grid phase 
spaces, shown in transparency, at uipet — 32. 
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Figure 9: Phase spaces after saturation for the Muhi Level Muhi Domain 
simulation of the Weibel instability, (a) Phase space Vx/c vs x/de (al) and Vy/c 
vs x/de (all) for the coarser grid and (b) phase space v^/c vs x/de (bl) and 
Vy/c vs x/de (bll) for the refined grid superimposed to the coarser grid phase 
spaces, shown in transparency, at uipet = 48. 
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Figure 10: Electric field mismatch between the levels in the overlapping area 
for the Multi Level Multi Domain simulation of the two stream instability, (a) 
Electric field E^.p^gig on the coarser grid alone and (b) with superimposed the 
refined grid field Ex^N,gn 3.96 — 0.5 < x/d^ < 6.95 + 0.5, with x/de = 3.96 and 
x/de = 6.95 being the refined grid boundaries. The contour of the coarser level 
density Ug^^ is superposed to both panels. The red rectangle and, in particular, 
the red line at ujpet = 62 mark an area of rather pronounced differences between 

Ex,P,gio and Ea;^N,gn- 
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Figure 11: Phase spaces at a time of rather pronounced differences between 
the longitudinal electric fields in the different levels for the Multi Level Multi 
Domain simulation of the two stream instability, (a) Phase space v^/c vs x/de 
(al) and Vy/cvs x/de (all) for the coarser grid and (b) phase space Vx/cvs xjdf, 
(bl) and Vyjc vs x/de (bll) for the refined grid superimposed to the coarser grid 
phase spaces, shown in transparency, at ujpet = 62. 



39 



(bl) 



0.5 







-0.5 



(bll) 



0.1 



-0.1 



Grid 0, (0 t= 92 

pe 



Superimposed, o>^t= 92 




10 



Superimposed, m t= 92 



2 4 6 

x/d 



10 




I 
I 



10 



-10 



Figure 12: Phase spaces at a time of minimal differences between the longi- 
tudinal electric fields in the different levels for the Multi Level Multi Domain 
simulation of the two stream instability, (a) Phase space v^/c vs x/d^ (al) and 
Vy/c VS x/de (all) for the coarser grid and (b) phase space Vx/c vs x/d^ (bl) 
and Vy/c vs x/de (bll) for the refined grid superimposed to the coarser grid 
phase spaces, shown in transparency, at uipet = 92. 
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Figure 13: Electric field Ey evolution when a shock wave is excited at the 
boundaries of the coarse domain. Panel (a) and (b) depict a reference simulation 
with the resolution of the coarser and refined grid respectively, panel (c) shows 
^v,N,gii superimposed to Ey^p^g^^ in a Multi Level Multi Domain simulation of 
shock propagation. 
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Figure 14: Current Jy evolution at 7 < a;/(ii < 13 when a shock wave is excited 
at the boundaries of the coarse domain. Panel (a) and (b) depict a reference 
simulation with the resolution of the coarser and refined grid respectively, panel 
(c) shows Jy,gii superimposed to Jy.giO in a Multi Level Multi Domain simulation 
of shock propagation. 
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